Nonadiabatic wavepacket dynamics: fc-space formulation 
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The time evolution of wavepackets in crystals in the presence of a homogeneous electric field 
is formulated in fc-space in a numerically tractable form. The dynamics is governed by separate 
equations for the motion of the waveform in fc-space and for the evolution of the underlying Bloch-like 
states. A one-dimensional tight-binding model is studied numerically, and both Bloch oscillations 
and Zener tunneling are observed. The long-lived Bloch oscillations of the wavepacket center under 
weak fields are accompanied by oscillations in its spatial spread. These are analyzed in terms of a 
fc-space expression for the spread having contributions from both the quantum metric and the Berry 
connection of the Bloch states. We find that when sizeable spread oscillations do occur, they are 
mostly due to the latter term. 

PACS numbers: 72.10.Bg, 78.20.-e, 78.20.Bh, 71.90.+q 



I. INTRODUCTION 

The study of the dynamics of electron wavepackets in 
crystals has experienced a revival in recent years. The 
development of heterostructure superlattices, photonic 
crystals, and optical lattices has opened new possibilities 
for the experimental realization of fundamental dynam- 
ical effects such as Bloch oscillations^^^ and Zener 
tunnclingi&i^ The wavepacket picture of transport has 
also shed light on subtle transport phenomena in solids. 
For instance, the intrinsic anomalous Hall effect in fer- 
romagnets was shown to result from a Berry-curvature 
term in the wavepacket group velocity^ 

Numerical simulations provide valuable insights into 
the dynamics of wavepackets in crystals. For instance, 
Bouchard and Lubani^ carried out a detailed study on 
a one-dimensional biased lattice, finding a rich variety 
of dynamical phenomena (Bloch oscillations of the cen- 
ter of mass, coherent breathing modes, Zener tunneling, 
and intrawell oscillations) as a function of the field-free 
band structure, field strength, and the form of the ini- 
tial wavepacket. In order to solve numerically the time- 
dependent Schrodingcr equation, they employed a super- 
ccll geometry with hard- wall boundary conditions; care 
had to be taken to ensure that the wavepacket never 
came close to the hard-wall boundaries for the duration 
of the simulation. In situations where unbounded accel- 
eration (via Zener tunneling) of a significant portion of 
the wavepacket takes place, a large supercell must then 
be used, which may become computationally demand- 
ing. In principle that can be avoided by switching from 
hard-wall to periodic boundary conditions. However, the 
inclusion in the Hamiltonian of the nonperiodic electric 
field term e£ ■ r then becomes problematic. A successful 
numerical strategy for describing homogeneous electric 
fields under periodic boundary conditions was developed 
in Refs. [Illfl2l for static fields, and generalized to time- 
dependent fields in Ref. [TH 

In Refs. [Tllfl2lfl3l the goal was to solve for the electronic 
structure of insulators in the presence of a homogeneous 



field. In this work we use a similar strategy to describe 
wavepacket dynamics. Our starting point is to express 
the wavepacket as a linear superposition of Bloch states, 

\4>)= [ ' dkf k \i/> k ) = f ' dke lkx f k \v k ) (1) 
Jo Jo 

(for simplicity we shall work in one dimension), and to 
follow the time evolution of both the waveform f k and the 
underlying states \v k )- If a wavepacket, initially prepared 
in a given band, is constrained to remain in the same 
band at later times, one obtains the "semiclassical" ap- 
proximation, which becomes exact in the adiabatic limit. 
Instead, we will allow for a fully unconstrained time evo- 
lution. As a result, for t > the states \ijjk) may become 
an admixture ^2 n c nk \ip^l) °f several eigenstates of the 
crystal Hamiltonian H°. We shall refer to such nonadi- 
abatic states as Bloch-like, since they retain the Bloch 
form, with V k (x + a) = v k {x). 

The manuscript is organized as follows. The equations 
of motion for f k and \v k ), and the expressions for the 
center and spread of the packet in terms of them, are 
derived in Sees. |H] and |TTT] respectively. In Sec. [TV] we 
perform simulations for a one-dimensional tight-binding 
Hamiltonian, using a numerically tractable form, given 
in the Appendix, of those equations. 



II. DYNAMICAL EQUATIONS IN fc-SPACE 

The wavepacket evolves according to the Schrodingcr 
equation (e = Ti = 1) 

4r = (*° ( 2 ) 

where £ is the electric field, which can be time-dependent 
but must be spatially uniform. Since fk and \vk) enter 
Eq. (UJ as a product, they are individually defined only 
up to a multiplicative factor which, because (vhlvh) = 1, 
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must take the form e lVk . We fix this phase arbitrariness 
by choosing f k to be real and positive. 
Inserting Eq. (fT]) into Eq. ©. 



i I dke lkx (f k \v k )+f k \v k )) = 

dke* kx H° k f k \v k )+£ [ dke tkx xf k \v k ), (3) 



where H k = e lkx fp e tkx an d henceforth the integration 
range from to 27r/a will be implied. Rewriting 

£ J dk(-id k e ikx )f k \v k ) 

= £ f dke lkx i (\v k )8 k f k + f k \d k v k )) (4) 



(d k = d/dk and an integration by parts was performed) 
then yields, at each k, 

ifk\v k ) +if k \v k ) = Hlf k \v k ) + 

+ i£(\v k )d k f k + f k \d k v k )). (5) 

Contracting with {v k | and subtracting from the resulting 
equation its complex conjugate, we arrive at the equation 
of motion for f k , 



fh — £d k f kl 



(0) 



where the reality of f k was used, together with the re- 
lations (v k \v k ) = -(v k \v k ) and (d k v k \v k ) = -(v k \d k v k ). 
To find the equation of motion for \v k ) we plug J6|) back 
into ©: 



i\v k ) = (H° k +i£d k ) \v k ). 



(7) 



Eqs. ©-([Jj govern the coherent wave pac ket dynamics. 
Eq. ([7]) was previously obtained in Ref. Il3l . where it was 
shown to describe the dynamics of valence electrons in 
insulators under the homogeneous field £(t). Here it de- 
scribes the nonadiabatic evolution of the Bloch-likc states 
\v k ) supporting the wavepacket. As for Eq. ©, it is the 
familiar result for the fc-space dynamics of the waveform, 
which remains valid in the presence of interband mixing. 

For numerical implementation the fc-derivatives must 
be replaced by finite-difference expressions over a fc-point 
mesh. While such discretization is straightforward for 
Eq. ©, Eq. ([7]) requires some care. As in Ref. [TJ|, we 
replace it with 



i\v' k ) = {H° k +i£d k )\v' k ), 



(8) 



where \d k v' k ) =Q k \d k v' k ) (Q k = 1 - \v k )(v k \ = 1 - P k ). 

Unlike \d k v k ), \d k v' k ) lends itself to a numerically robust 
finite-differences representation (see Appendix). 

The states \v' k ) obeying Eq. ([8]) differ from the states 
\v k ) in Eq. (UJ by a phase factor, 



which we must keep track of. Inserting ^ into © and 
using ([Jj yields 



iU k = -£A k U k , 
where A k is the Berry connetion, 

A k = i(v k \d k v k ). 



(10) 



(11) 



Eqs. ©, and (flTj)) are the desired dynamical equa- 
tions for f k , \v' k ), and U k . Together they determine the 
time evolution of the wavepacket 



dkf k e ikx U* k \v' k ). 



(12) 



In practice Eqs. ([8]) and (fTTj) arc replaced by the dis- 
cretized forms (|A.3[) and (|A.6|) . respectively. 



III. WAVEPACKET CENTER AND SPREAD 

In the previous Section we formulated the dynamics 
of the wavepacket (UJ in terms of f k and \v k ). Here we 
shall express its center and spread in terms of those same 
fc-spacc quantities. 



A. fc-space expressions 

Let us define the generating function for the spatial 
distribution of the wavepacket: 

C(q) = = — / dkf k f k+q (v k \v k+q ), (13) 



where the second equality follows from Eq. (Q} together 
with the identity 



{^ kl |e"^|^ 2 ) = — S(ka-ki-q)(v hl \ Vkl+q ). (14) 
The first moment is given by 

(x) = id q C(q)\ q=Q = (%D h ) = {A k ), (15) 
where iD k is the Hermitian operator 

iD k =id k + A k (16) 
and we have introduced the notation 



2-7T f 

{O k ) = —J dkf k O k f k 



(17) 



The last equality in Eq. (jT5]) follows from f k being real. 
Next we evaluate the spread 

(Axf = (x 2 ) - {x)\ (18) 

For (x) 2 we use (fl~5|) . while (x 2 ) is given by i 2 d q C(q)\ q= o- 

2tt 



Vk) = U k \v k ), 



(9) 



(' ) = 



dk(d k / k y+ / dk f£(d k v k \d k v k ) 



(19) 
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Inserting 1 = P k + Q k in the last term on the RHS and 
then combining with Eq. (fT5|) yields 

(Ax) 2 = ^J dk(8 k f k ) 2 + (G k ) + {(AA k ) 2 ) , (20) 

where G k is the quantum metric^ 

G k = (d k v k \d k v k ). (21) 

All three terms in Eq. (|20|) arc manifestly non-negative. 
The first one only depends on f k , while the remaining two 
also depend on the states \v k ). However, the second term 
is insensitive to the phases of those states (it is invariant 
under the gauge transformation (|A.9|) ). whereas the third 
term is phase-dependent. 

It is instructive to consider the limit of a uniform 
waveform, f k = a/2ir, in which \<f>) becomes a Wannier 
function. Eq. (fT5")) can then be recast as (x) = cnp/2ir, 
where <p = J dk A k is the Berry phase associated with 
the manifold of states \v k )J£- The first term on the RHS 
of Eq. ([20]) then vanishes identically, while the second 
and third terms reduce to the gauge-invariant and gauge- 
dependent parts of the Wannier spread for an isolated 
band in one dimension, given respectively by Eqs. (C12) 
and (C17) of Ref. El 

B. Uncertainty relation and minimal wavepackets 

An alternative decomposition of the wavepacket spread 
may be obtained by noting that 

( (tD k ) 2 ) = ^ J dk(d k f k ) 2 + (A 2 ), (22) 

as can be readily verified using the hermiticity of iD k 
and the reality of f k . Comparison with Eqs. (fl~5|) and 
(|2H|) shows that 

(Ax) 2 = {(A(iD k )) 2 ) + (G k ). (23) 

Combining this with the relation 

(AA) 2 (AB) 2 >l\{[A,B})\ 2 (24) 

yields, upon setting A = iD k , B = k, and using 

[iD k ,k] = i, 

[(Ax) 2 -(G k )](Ak) 2 >± (25) 

In the limit of a vanishing lattice potential G k — > and 
k — > p (canonical momentum). Eq. ([25]) then reduces the 
familiar Heiscnbcrg uncertainty relation. 

Let us now show that (|25[) becomes an equality for 
minimal wavepackets in one dimension. Once the mani- 
fold of states \v k ) and the width Afc of the waveform are 
specified, all that remains is to set the phases of the |ufc) 
and the shape of f k . We wish to minimize the spread (|2"0|) 



with respect to those two parameters. We start with the 
phases, which only affect the term ((AA k ) 2 ). This term 
vanishes when A k is constants in which case 

2-7T f 

(Ax) 2 - (G k ) = — I dk(d k f k ) 2 . (26) 
Combining the previous two equations, 

~a J dk ^ h)2 > < 27 > 

It can be verified that for Afc <C 2ir/a this becomes an 
equality when the waveform has a Gaussian shape. In 
conclusion, a minimal wavepacket in one dimension is 
characterized by a Gaussian-shaped f k and a constant 
Berry connection A k in the region of fc-space where f k is 
non-negligible. Its spread equals 

(Ax) 2 lin = 47^2 + < G *>- ( 28 ) 

It was shown in Ref. that the spread of a maximally- 
localized Wannier function in one dimension is (Ax)^ in = 
(Gk)- This result can be viewed as the limit Afc — * 00 
of Eq. (|28p . In the opposite limit of a narrow waveform, 
the wavepacket spread becomes dominated by the term 
l/4(Afc) 2 , as will be illustrated in the next Section. 

IV. NUMERICAL RESULTS 

A. Tight-binding model 

We have applied our scheme to the same one- 
dimensional tight-binding model used in Ref. 115 . This 
is a three-band Hamiltonian with three atoms per unit 
cell of length a = 1 and one orbital per atom, 

H ° = E {'.*<v:> + 7 [*h+i + 4+1*} } . ( 29 ) 

3 

with the site energy given by €3™+; = U cos /3; . Here 
m is the cell index, I = { — 1,0, 1} is the site index, and 
0i = 2ttI/3. The upper panel of Fig. Q] shows the energy 
dispersion for 7 = — U =1. 

Before the spatial distribution of the wavepacket can 
be defined, the matrix elements of the position operator 
must be specified. As in Ref. [H we choose the simplest 
diagonal representation x = Ylj x j c j c j> with Xj = j/3. 
The lower panel of Fig. [1] shows the quantum metric cal- 
culated for the Bloch states in the lowest band. As ex- 
pected from the relation^ G k < 1/2 A k (A k is the direct 
gap to the second band), G k peaks around k = 0. 

We will now study numerically the wavepacket dy- 
namics on this model. In Sees. IIVBI and IIV CI we con- 
sider respectively Bloch oscillations in a weak field and 
Zener tunneling in a strong field. In both cases we 
gradually turn on the field over a time interval, T, as 
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FIG. 1: Upper panel: Band structure of the one-dimensisonal 
tight-binding model of Eq. (|29p . for the choice of parameters 
7 = — J7 = 1. Lower panel: Quantum metric [Eq. (|21|) ] for 
the lowest band. 



£(t) = £o sin(7rf/2T), and keep it constant afterwards. 
The simulations begin with a minimal wavepacket pre- 
pared in the lowest band. Unless otherwise noted, the 
width of the Gaussian waveform is Afc = 0.075 x 2tt. 



B. Bloch oscillations in a weak electric field 

In order to observe long-lived Bloch oscillations we 
choose a weak field £q = 0.055 Ao/o, which we turn 
on over a time interval of the order of the Bloch oscilla- 
tion period tb = 2Tr/£oa (henceforth in this subsection 
we choose t = long after the field has saturated at So). 
100 fc-points are used to sample the Brillouin zone, and 
the time step is At = 1.7 x 10~ 5 tb. 

In the upper panel of Fig. [5] the weights \{j\4>)\ 2 of 
the wavepacket on the tight-binding orbitals \j) are used 
to depict its spatial distribution as a function of time. 
The Bloch oscillations of (x) are clearly seen. In the 
lower panel we plot the oscillation amplitude, A, ver- 
sus the width, W, of the lowest band, which was tuned 
by adjusting the tight-binding parameters in the range 
0.5 < 7 < 1.5 and -1.5 < U < -0.5. A and W are lin- 
early related, with a slope A/W ~ 8.5. This is in agree- 
ment with the prediction A = (1/2)(W/£o)\Sq\, valid for 
a single-band tight-binding modelJ^ The dimensionless 
parameter Sq depends on the initial choice of wavepacket, 
but its magnitude cannot exceed 1; in the present case 
l/2£ - 9.09, so that \S \ ~ 0.9. 

Next we analyze the behavior of the wavepacket spread 
in the course of the Bloch oscillations, by tracking each of 
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FIG. 2: (Color online.) Bloch oscillations of a wavepacket 
prepared in the lowest band. Upper panel: Time evolution of 
|(j|0)| 2 , the weights of the wavepacket on the tight-binding 
basis orbitals, for tight-binding parameters 7 = — U = 1. 
Lower panel: Amplitude A of the Bloch oscillations versus 
the bandwidth W. 



the three terms in Eq. (f20|) . According to Eq. ([6|), as the 
center (k) of the packet traverses the Brillouin zone, the 
shape of the waveform fk remains unchanged. Hence the 
term (2-k/o) J dk(dkfk) 2 is a constant of motion. Oscil- 
lations in the spread must therefore arise from the other 
two terms, (G) and ((AA) 2 ) (henceforth the k subscript 
will be ommited for brevity). 

Since initially the wavepacket was minimal 
(((AA) 2 ) = 0), one might have expected spread 
oscillations to arise mostly from the fc-space dispersion 
of the metric, i.e., from the varying constraint imposed 
by the uncertainty relation (|25p as the wavepacket moves 
through fc-space. Instead we find that after an initial 
transient the spread of the packet can become far from 
minimal, with ((A^4) 2 } 3> (G) (the bars denote a time 
average over several Bloch oscillations). This is the 
situation depicted in the upper panel of Fig. [31 which 
pertains to the same simulation run as in the upper 
panel of Fig. [21 Note that the spread undergoes sizeable 
oscillations, arising mostly from ((A^4) 2 ) (not shown). 

The above scenario is typical of Bloch oscillations in a 
relatively wide band. When we adjust the tight-binding 
parameters so as to reduce the bandwidth, we find that 
the wavepacket remains close to minimal in the course 
of the Bloch oscillations, and the associated spread oscil- 
lations are very weak, of the order of the Brillouin zone 
dispersion of the metric (lower panel of Fig.[3J). One could 
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FIG. 4: Dependence of the time-averaged contributions to the 
wavepacket spread [Eq. (|20p ] on the waveform width Afc, for 
7 = — U = 1. The two dashed lines are fits to the values 
from the numerical simulation (symbols), while the solid line 
equals l/4(Afc) 2 , the right-hand-side of Eq. ((2"T)l . 



FIG. 3: Time evolution of the total wavepacket spread (Ax) 2 
(solid lines), and of the metric contribution {Gk} (dashed 
lines), for two different sets of tight-binding parameters. Up- 
per panel: 7 = — U = 1, resulting in a bandwidth W ~ 0.5 
and a minimum bandgap Ao ~ 1.1. Lower panel: 7 = 0.5 
and U = -1, for which W ~ 0.1 and A ~ 1.2. 



try to enhance the spread oscillations by further tuning 
the model parameters so as to make the metric very large 
in some regions of fc-space. However, that would require 
very small gaps, and under those circumstances the Bloch 
oscillations are strongly damped by Zener tunneling. 

So far we have considered a fixed waveform width 
Afc = 0.075 x 2ir. Fig. 0] displays the dependence on 
Afc of each contributions to (Ax) 2 , averaged over sev- 
eral Bloch oscillations. The term (G) is roughly constant 
and equal to the Brillouin zone average of the metric; it 
remains a small fraction of (Ax) 2 over the entire range 
of Afc. For Afc <C 2ir/a the spread is dominated by the 
term (2-k / a) J dk(dkfk) 2 , while for larger values of Afc 
the term ((AA) 2 ) takes over. Its monotonic increase is 
easily understood: the larger the range Afc, the larger 
the spread of Ak over that range is likely to be. 



C. Zener tunneling in a strong electric field 




FIG. 5: (Color online.) Upper panel: Time evolution of a 
wavepacket prepared in the lowest band, with child packets 
in the second and third bands appearing as a result of Zener 
tunneling. The tight-binding parameters are 7 = 1 and U = 
—0.3. Lower panel: Electric field as a function of time. The 
field is switched on from t = to t = 0.86tb, where tb is the 
Bloch period for the saturation field £q. 



In the previous subsection a weak electric field (£q 
Ao /a) was chosen, so that for the duration of the simula- 
tion the wavepacket remained mostly in the lowest band. 
For sufficiently strong fields, significant interband transi- 
tions are expected to occur as the wavepacket reaches the 
zone center, where the gap is smallest. In order to ob- 
serve this phenomenon! the saturation field was increased 



from Eq = 0.055 to £q = 0.09, while the minimum gap 
between the first two bands was reduced from Ao = 1.14 
to A = 0.31 by setting U = -0.3. With this choice of 
parameters a very good stability of the propagation al- 
gorithm is needed, and we decreased the time step from 
2 x 10~ 3 to 2 x 10~ 4 . The Brillouin zone was sampled 
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FIG. 6: Band occupancy P n (t) for the simulation shown in 
Fig. [5] 



over 300 £;-points. 

The Zener tunneling can be seen in the upper panel of 
Fig.[5]as a splitting of the wavepacket in real space, At the 
end of every Bloch oscillation, t = n ■ tb (n = 1, 2, . . . ), 
the main wavepacket in the lowest band spawns child 
packets which oscillate in the second band with the same 
period tb but a larger amplitude (due to the larger width 
of the second band) and gain more weight after each 
Bloch cycle. These child packets also spawn grandchild 
packets at t = (2n + 1)tb/2, when Zener tunneling from 
the second to the third band becomes possible. 

A more quantitative picture is obtained by monitoring 
the distribution of the packet among the three bands. We 
define the band occupancy P n (t) as the total probability 
that the wavepacket resides on band n-J£ 



Pn{t) = 



E 

k 



fl{t)\{u ( % k (t))\ 



(30) 



V. CONCLUSIONS 

In this work we have developed a new numerical 
scheme for simulating wavepacket dynamics in a peri- 
odic lattice potential with a linear potential (homoge- 
neous electric field) superimposed. By using the A:-space 
representation of the position operator, we were able to 
include the linear electric-field term e£ ■ r in the Hamilto- 
nian under periodic boundary conditions, thus avoiding 
having to use large supercells with hard-wall boundary 
conditions^ 

In the present approach the wavepacket is represented 
on a uniform mesh of fc-points by a waveform sitting 
on top of a "band" of states \vk) [Eq. ([T])]. The time 
evolution of the wavepacket is then obtained from that 
of fk and \vk). For £ ^ the states \vk) become non- 
adiabatic, field-polarized Bloch states which span several 
energy bands of the field-free crystal Hamiltonian H° (a 
similar representation of wavepackets in coupled energy 
bands was used in Ref. [l8h : thus interband effects such 
as Zener tunneling are fully accounted for. 

The method was tested on a one-dimensional tight- 
binding model. Depending on the choice of tight-binding 
parameters and electric field strength, we observed cither 
long-lived Bloch oscillations, or short-lived Bloch oscil- 
lations strongly damped by Zener tunning. In the for- 
mer regime we monitored the changes in the wavepacket 
spread accompanying the Bloch oscillatations of the cen- 
ter of mass, and identified two distinct situations: (i) For 
wavepackets moving in narrow bands, the spread changed 
very little over time, (ii) For wavepackets moving in wide 
bands, the Bloch oscillations were accompanied by con- 
siderable oscillations of the wavepacket spread. An anal- 
ysis of the fc-space expression for the spread [Eq. (|20[) ] 
reveals two distinct contributions which can change over 
time: one associated with the Berry connection of the un- 
derlying Bloch states, and another related to the quan- 
tum metric. By tracking each of them separately, we 
concluded that in the cases where significant spread os- 
cillations took place, they originated mostly in the Berry 
connection term. 



This quantity is plotted in Fig. [6] as a function of time. 
Initially only the first band is occupied. After a Bloch 
period, the wavepacket reaches the zone center, where 
the gap to the second band is smallest, at which point 
significant Zener tunneling occurs, giving rise to a partial 
occupation of the second band. (Between t = and t = 
tb the wavepacket moved in fc-space by less than the full 
Brillouin zone width, from (k) = —2tt/3 to (k) ~ — 2ir, 
because during most of that time interval the electric field 
strength was less than as shown in the lower panel 
of Fig. ) Subsequently there are also transitions to the 
third band, again with periodicity tb] because transitions 
from the second band to the first and third bands happen 
at the zone center and at the zone boundary respectively, 
Pi undergoes changes twice as often as P\ and P3. 
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APPENDIX: DISCRETIZED EXPRESSIONS IN 
fc-SPACE 

Here we derive the discrctized versions used in Sec. IIVI 
of the dynamical equations of Sec. |TT] and of the 
wavepacket center and spread expressions of Sec. IIIII 
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1. Dynamical equations 



use Eq. (|A.6[) . A few manipulations yield 



The appropriate finite-difference representation of 
\dkv' k ) on a uniform grid is^ 3 - 



\dkV k ) = ^(\ v k+b) - \Vk-b}): 
where b is the mesh spacing and 

K +6 > 



\Vk+b) 



We use Eq. ([AT]) to recast Eq. © a 
i\v' k )=T k \v' k ), 
in terms of the Hcrmitian operator 



13 



(A.l) 



(A.2) 



(A.3) 



T k = Hi 



w 



fcj 



where w k = (i£/2b) (P+ - P fe ") and P fe d 



(A.4) 

l«fe±6)KI- 

Eq. (|A.3|) is solved numerically at each grid point 

i(At/2)T fc (t). 



usin g 10 ' 13 



KX* + At) 



l + i(At/2)T fc (t) 



!<)(*)• 



(A.5) 



Because of the similarity between Eqs. (|T0|) and (|A.3|) . 
the phase factors £/fc can be propagated in time using the 
same algorithm. A finite-difference representation for the 
connection in Eq. (fT0|) is needed. We usei^ 



2b V fe 



$7 



(A.6) 



where we have defined the phases $^ = — lmhi{v k \v k ±b). 
When propagating U k via Eq. (fTU|) , care must be taken in 
choosing consistently the branch cuts for the two phases 
in Eq. (|A.6[) . to ensure that A k remains a smooth function 
of k at every time step. 



2n 



a 2 

k 



b) $ fc ■ 



(A.7) 



To find an expression for (Ax) 2 we start from Eq. (|20[) . 
The first term on the RHS is easily discretized, and (G k ) 
can be evaluated from (f2"Tj) . The remaining term is (Ai) — 
(A k ) 2 . For (A k ) 2 = (x) 2 we use (|A.7|) and, from (|A.6j) . 



/fc + fk+b , 



(A.8) 



Besides reducing to the correct continuum expressions 
as b — > 0, Eqs. (|A.7[) and (|A.8[) preserve exactly, for finite 
&, certain properties of those expressions. If we perform 
a change of phases 



(A.9) 



with 9 k ~ "f k — kR (R is a lattice vector and 7fc+2u-/a 
7fc), the center of the packet, Eq. (fT5|) . changes as 



(a?) ^ (x) + R / dk f 2 d klk . (A. 10) 

a 



For a Wannier wavepacket is constant, so that the 
last term vanishes, and (x) changes at most by a lat- 
tice vector jiii If instead f k spans a narrow region /C of 
the Brillouin zone, (x) can then shift continuously under 
()A.9[) . Consider the choice 



-5k keIC 







otherwise ' 



(A.11) 



2. Wavepacket center and spread 

In order to obtain a finite-difference representation of 
Eq. (fT5)) we make the replacement J dk — > b and then 



which produces a rigid shift (x) — ■> (x) + (5. Eq. (|A. 7|) 
obeys this transformation exactly, while the spread eval- 
uated using (|A.8|) remains unchanged, as can be easily 
verified. 
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